library('tidyverse')
library('tidyquant')
library('tseries')
library('gghighlight')
load('api_key.rdata')
Time series analysis:
- Stationary vs. Non-Stationary
- Additive vs. Exponential
- ARMA / ARIMA models
- Filters: Kalman, Holt-Winters
- Stochastic Processes, Random Walk, Wiener Process, Brownian Motion
What is a stationary time series?
“A time series is stationary if a shift in time doesn’t cause a change in the shape of the distribution. Basic properties of the distribution like the mean, variance, and covariance are constant over time.”
Different types of stationarity:
- Strict Stationarity - basically doesn’t exist… essentially nothing changes based off of time…
- First-Order Stationarity - means never change, but other statistics can change
- Second-Order (Weak) Stationarity - constant mean, variance, and autocovariance that doesn’t change with time
- Trend-Stationary Models - fluctuate around deterministic trend (series mean). Trend is typically linear or quadratic, amplitude of fluctuations stays constant
- Difference-Stationary Models - one or more differences of the data need to be taken to make it stationary
Why do we care about stationarity?
Most statistical models that have been developed assume stationarity.
set.seed(100)
# monthly data, normally distributed
# 100 stoltzcoin average trade price per month
# Standard deviation of 3
dates = seq(from = lubridate::ymd('2011-07-01'),
to = lubridate::ymd('2018-06-30'),
by = 'months')
n = length(dates)
mu = 100
sd = 3
time_series = tibble(date = dates, stoltzcoin = rnorm(n = n, mean = mu, sd = sd),
type = 'stationary', trend = 'none')
ggplot(time_series, aes(x = date, y = stoltzcoin)) +
geom_line() + geom_smooth(method = lm,
formula = y ~ x, se = FALSE) +
geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Always test for stationarity even if it looks stationary
Augmented Dickey-Fuller is a typical test - tseries::adf.test()
- Tests null hypothesis that the autoregressive model has a unit root
- Stationarity means rejecting the null hypothesis
- We will look for p-value < 0.05
# adf.test only accepts univariate time series data, zoo object will be used
adf.test(zoo::zoo(time_series$stoltzcoin, time_series$date))
p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: zoo::zoo(time_series$stoltzcoin, time_series$date)
Dickey-Fuller = -4.1967, Lag order = 4, p-value = 0.01
alternative hypothesis: stationary
Add a linear trend for the example
Stoltzcoin has been going up by 2 each month
Known as “trend stationary”
going_up_by = 2
trend = seq(from = 0,
to = (n-1) * going_up_by,
by = going_up_by)
ts_lin_trend = tibble(date = dates, stoltzcoin = time_series$stoltzcoin + trend,
type = 'non-stationary', trend = 'linear')
ggplot(ts_lin_trend, aes(x = date, y = stoltzcoin)) +
geom_line() + geom_smooth(method = lm,
formula = y ~ x, se = FALSE) +
geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Add a seasonal component to the linear trend
Stoltzcoin fluctuates throughout the year
going_up_by = 2
trend = seq(from = 0,
to = (n-1) * going_up_by,
by = going_up_by)
seasonal = sin(1:n) * 10
ts_lin_trend_seasonal = tibble(date = dates, stoltzcoin = time_series$stoltzcoin + trend + seasonal,
type = 'non-stationary', trend = 'linear')
ggplot(ts_lin_trend_seasonal, aes(x = date, y = stoltzcoin)) +
geom_line() + geom_smooth(method = lm,
formula = y ~ x, se = FALSE) +
geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

Add an exponential trend to the seasonal component
Stoltzcoin has increased in volatility over the years
- Stoltzcoin now has:
- A positive and growing linear trend
- Seasonality
going_up_by = 2
trend = seq(from = 0,
to = (n-1) * going_up_by,
by = going_up_by)
seasonal = sin(1:n) * 1:n
ts_lin_trend_seasonal_exp = tibble(date = dates, stoltzcoin = time_series$stoltzcoin + trend + seasonal,
type = 'non-stationary', trend = 'linear')
ggplot(ts_lin_trend_seasonal_exp, aes(x = date, y = stoltzcoin)) +
geom_line() + geom_smooth(method = lm,
formula = y ~ x, se = FALSE) +
geom_hline(aes(yintercept = mean(stoltzcoin)), color="red")

R’s decompose function does a lot for you
d1y = as.numeric(lubridate::year(min(ts_lin_trend_seasonal_exp$date)))
d1m = as.numeric(lubridate::month(min(ts_lin_trend_seasonal_exp$date)))
d2y = as.numeric(lubridate::year(max(ts_lin_trend_seasonal_exp$date)))
d2m = as.numeric(lubridate::month(max(ts_lin_trend_seasonal_exp$date)))
x = ts(ts_lin_trend_seasonal_exp$stoltzcoin, start=c(d1y, d1m), end=c(d2y, d2m), frequency=12)
plot(decompose(x, type = 'additive'))

plot(decompose(x, type = 'multiplicative'))

Example from decompose function
random portion should just look like white noise (otherwise you’ve likely missed something)
require(graphics)
m <- decompose(co2)
plot(m)

dow = tq_index("DOW") %>%
tq_get(get = "stock.prices")
head(dow, 10)
p = ggplot(dow, aes(x = date, y = close, col = symbol)) +
geom_line() +
gghighlight(max(close) > 250) +
theme_minimal() +
facet_wrap(~ symbol)
print(p)

Is BA a stationary time series?
Can you simply zoom in and find a portion that is stationary?
ba = dow %>%
filter(symbol == 'BA') %>%
select(date, close)
plotly::ggplotly(ggplot(ba, aes(x = date, y = close)) + geom_line() + geom_smooth(method = lm, se = FALSE))
ts_fit
Call:
lm(formula = close ~ date, data = ba)
Coefficients:
(Intercept) date
-734.75085 0.05385


Are the daily returns stationary?
ret = ba %>%
tq_transmute(close, mutate_fun = dailyReturn)
ggplot(ret, aes(x = date, y = daily.returns)) + geom_line() + geom_smooth(method = lm, se = FALSE)

adf.test(zoo::zoo(ret$daily.returns, ret$date))
p-value smaller than printed p-value
Augmented Dickey-Fuller Test
data: zoo::zoo(ret$daily.returns, ret$date)
Dickey-Fuller = -14.308, Lag order = 13, p-value = 0.01
alternative hypothesis: stationary
ts_fit = lm(daily.returns ~ date, data = ret)
ts_fit
Call:
lm(formula = daily.returns ~ date, data = ret)
Coefficients:
(Intercept) date
-8.376e-03 5.736e-07
res = tibble(residuals = resid(ts_fit))
ggplot(res, aes(x = 1:nrow(res), y = residuals)) + geom_point() + geom_rug(sides = 'l')

ggplot(res, aes(x = residuals)) + geom_histogram(aes(y = ..density..), bins = 40) + geom_density(col = 'red', size = 1.5)

LS0tCnRpdGxlOiAiVGltZSBTZXJpZXMgTW9kZWxpbmcsIEFuYWx5c2lzIGFuZCBEZWNvbXBvc2l0aW9uIgpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sKLS0tCgpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0KbGlicmFyeSgndGlkeXZlcnNlJykKbGlicmFyeSgndGlkeXF1YW50JykKbGlicmFyeSgndHNlcmllcycpCmxpYnJhcnkoJ2dnaGlnaGxpZ2h0JykKbG9hZCgnYXBpX2tleS5yZGF0YScpCmBgYAoKCiMjIyMgVGltZSBzZXJpZXMgYW5hbHlzaXM6ICAKCiAgKiBTdGF0aW9uYXJ5IHZzLiBOb24tU3RhdGlvbmFyeQogICogQWRkaXRpdmUgdnMuIEV4cG9uZW50aWFsCiAgKiBBUk1BIC8gQVJJTUEgbW9kZWxzCiAgKiBGaWx0ZXJzOiBLYWxtYW4sIEhvbHQtV2ludGVycwogICogU3RvY2hhc3RpYyBQcm9jZXNzZXMsIFJhbmRvbSBXYWxrLCBXaWVuZXIgUHJvY2VzcywgQnJvd25pYW4gTW90aW9uCiAgCgojIyMjIFdoYXQgaXMgYSBzdGF0aW9uYXJ5IHRpbWUgc2VyaWVzPwoKIkEgdGltZSBzZXJpZXMgaXMgc3RhdGlvbmFyeSBpZiBhIHNoaWZ0IGluIHRpbWUgZG9lc24ndCBjYXVzZSBhIGNoYW5nZSBpbiB0aGUgc2hhcGUgb2YgdGhlIGRpc3RyaWJ1dGlvbi4gQmFzaWMgcHJvcGVydGllcyBvZiB0aGUgZGlzdHJpYnV0aW9uIGxpa2UgdGhlIG1lYW4sIHZhcmlhbmNlLCBhbmQgY292YXJpYW5jZSBhcmUgY29uc3RhbnQgb3ZlciB0aW1lLiIKCiFbU3RhdGlvbmFyeSB2cyBTZWFzb25hbCBUaW1lIHNlcmllcyAtIGh0dHA6Ly93d3cuc3RhdGlzdGljc2hvd3RvLmNvbS9zdGF0aW9uYXJpdHkvXShpbWFnZTEucG5nKQoKCgojIyMjIERpZmZlcmVudCB0eXBlcyBvZiBzdGF0aW9uYXJpdHk6ICAKCiAgKiBTdHJpY3QgU3RhdGlvbmFyaXR5IC0gYmFzaWNhbGx5IGRvZXNuJ3QgZXhpc3QuLi4gZXNzZW50aWFsbHkgbm90aGluZyBjaGFuZ2VzIGJhc2VkIG9mZiBvZiB0aW1lLi4uCiAgKiBGaXJzdC1PcmRlciBTdGF0aW9uYXJpdHkgLSBtZWFucyBuZXZlciBjaGFuZ2UsIGJ1dCBvdGhlciBzdGF0aXN0aWNzIGNhbiBjaGFuZ2UKICAqIFNlY29uZC1PcmRlciAoV2VhaykgU3RhdGlvbmFyaXR5IC0gY29uc3RhbnQgbWVhbiwgdmFyaWFuY2UsIGFuZCBhdXRvY292YXJpYW5jZSB0aGF0IGRvZXNuJ3QgY2hhbmdlIHdpdGggdGltZSAKICAqIFRyZW5kLVN0YXRpb25hcnkgTW9kZWxzIC0gZmx1Y3R1YXRlIGFyb3VuZCBkZXRlcm1pbmlzdGljIHRyZW5kIChzZXJpZXMgbWVhbikuIFRyZW5kIGlzIHR5cGljYWxseSBsaW5lYXIgb3IgcXVhZHJhdGljLCBhbXBsaXR1ZGUgb2YgZmx1Y3R1YXRpb25zIHN0YXlzIGNvbnN0YW50CiAgKiBEaWZmZXJlbmNlLVN0YXRpb25hcnkgTW9kZWxzIC0gb25lIG9yIG1vcmUgZGlmZmVyZW5jZXMgb2YgdGhlIGRhdGEgbmVlZCB0byBiZSB0YWtlbiB0byBtYWtlIGl0IHN0YXRpb25hcnkKCgojIyMjIFdoeSBkbyB3ZSBjYXJlIGFib3V0IHN0YXRpb25hcml0eT8gIAoKTW9zdCBzdGF0aXN0aWNhbCBtb2RlbHMgdGhhdCBoYXZlIGJlZW4gZGV2ZWxvcGVkIGFzc3VtZSBzdGF0aW9uYXJpdHkuCgoKICAKYGBge3IsIHdhcm5pbmc9RkFMU0V9CnNldC5zZWVkKDEwMCkKCiMgbW9udGhseSBkYXRhLCBub3JtYWxseSBkaXN0cmlidXRlZAojIDEwMCBzdG9sdHpjb2luIGF2ZXJhZ2UgdHJhZGUgcHJpY2UgcGVyIG1vbnRoCiMgU3RhbmRhcmQgZGV2aWF0aW9uIG9mIDMKZGF0ZXMgPSBzZXEoZnJvbSA9IGx1YnJpZGF0ZTo6eW1kKCcyMDExLTA3LTAxJyksIAogICAgICAgICAgICB0byA9IGx1YnJpZGF0ZTo6eW1kKCcyMDE4LTA2LTMwJyksIAogICAgICAgICAgICBieSA9ICdtb250aHMnKQoKbiA9IGxlbmd0aChkYXRlcykKbXUgPSAxMDAKc2QgPSAzCgp0aW1lX3NlcmllcyA9IHRpYmJsZShkYXRlID0gZGF0ZXMsIHN0b2x0emNvaW4gPSBybm9ybShuID0gbiwgbWVhbiA9IG11LCBzZCA9IHNkKSwgCiAgICAgICAgICAgICAgICAgICAgIHR5cGUgPSAnc3RhdGlvbmFyeScsIHRyZW5kID0gJ25vbmUnKQoKZ2dwbG90KHRpbWVfc2VyaWVzLCBhZXMoeCA9IGRhdGUsIHkgPSBzdG9sdHpjb2luKSkgKyAKICBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgc2UgPSBGQUxTRSkgKyAKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gbWVhbihzdG9sdHpjb2luKSksIGNvbG9yPSJyZWQiKQpgYGAKCgoKIyMjIyBQZXJmb3JtIHJlZ3Jlc3Npb24gYW5kIGxvb2sgYXQgdGhlIHJlc2lkdWFscwoKYGBge3J9CnRzX2ZpdCA9IGxtKHN0b2x0emNvaW4gfiBkYXRlLCBkYXRhID0gdGltZV9zZXJpZXMpCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCiMjIEFsd2F5cyB0ZXN0IGZvciBzdGF0aW9uYXJpdHkgZXZlbiBpZiBpdCBsb29rcyBzdGF0aW9uYXJ5CiMjIyMjIEF1Z21lbnRlZCBEaWNrZXktRnVsbGVyIGlzIGEgdHlwaWNhbCB0ZXN0IC0gdHNlcmllczo6YWRmLnRlc3QoKQogICogVGVzdHMgbnVsbCBoeXBvdGhlc2lzIHRoYXQgdGhlIGF1dG9yZWdyZXNzaXZlIG1vZGVsIGhhcyBhIHVuaXQgcm9vdAogICogU3RhdGlvbmFyaXR5IG1lYW5zIHJlamVjdGluZyB0aGUgbnVsbCBoeXBvdGhlc2lzCiAgKiBXZSB3aWxsIGxvb2sgZm9yIHAtdmFsdWUgPCAwLjA1CgoKYGBge3J9CiMgYWRmLnRlc3Qgb25seSBhY2NlcHRzIHVuaXZhcmlhdGUgdGltZSBzZXJpZXMgZGF0YSwgem9vIG9iamVjdCB3aWxsIGJlIHVzZWQKYWRmLnRlc3Qoem9vOjp6b28odGltZV9zZXJpZXMkc3RvbHR6Y29pbiwgdGltZV9zZXJpZXMkZGF0ZSkpCgpgYGAKCgoKIyMjIEFkZCBhIGxpbmVhciB0cmVuZCBmb3IgdGhlIGV4YW1wbGUKIyMjIyBTdG9sdHpjb2luIGhhcyBiZWVuIGdvaW5nIHVwIGJ5IDIgZWFjaCBtb250aAojIyMjIEtub3duIGFzICJ0cmVuZCBzdGF0aW9uYXJ5IgoKYGBge3J9CmdvaW5nX3VwX2J5ID0gMgp0cmVuZCA9IHNlcShmcm9tID0gMCwgCiAgICAgICAgICAgIHRvID0gKG4tMSkgKiBnb2luZ191cF9ieSwKICAgICAgICAgICAgYnkgPSBnb2luZ191cF9ieSkKCnRzX2xpbl90cmVuZCA9IHRpYmJsZShkYXRlID0gZGF0ZXMsIHN0b2x0emNvaW4gPSB0aW1lX3NlcmllcyRzdG9sdHpjb2luICsgdHJlbmQsIAogICAgICAgICAgICAgICAgICAgICB0eXBlID0gJ25vbi1zdGF0aW9uYXJ5JywgdHJlbmQgPSAnbGluZWFyJykKCmdncGxvdCh0c19saW5fdHJlbmQsIGFlcyh4ID0gZGF0ZSwgeSA9IHN0b2x0emNvaW4pKSArIAogIGdlb21fbGluZSgpICsgZ2VvbV9zbW9vdGgobWV0aG9kID0gbG0sIAogICAgICAgICAgICAgICAgICAgICAgICAgICAgZm9ybXVsYSA9IHkgfiB4LCBzZSA9IEZBTFNFKSArIAogIGdlb21faGxpbmUoYWVzKHlpbnRlcmNlcHQgPSBtZWFuKHN0b2x0emNvaW4pKSwgY29sb3I9InJlZCIpCmBgYAoKCgojIyMjIFBlcmZvcm0gcmVncmVzc2lvbiBhbmQgbG9vayBhdCB0aGUgcmVzaWR1YWxzCgpgYGB7cn0KdHNfZml0ID0gbG0oc3RvbHR6Y29pbiB+IGRhdGUsIGRhdGEgPSB0c19saW5fdHJlbmQpCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCmBgYHtyfQphZGYudGVzdCh6b286Onpvbyh0c19saW5fdHJlbmQkc3RvbHR6Y29pbiwgdHNfbGluX3RyZW5kJGRhdGUpKQpgYGAKCgoKIyMjIEFkZCBhIHNlYXNvbmFsIGNvbXBvbmVudCB0byB0aGUgbGluZWFyIHRyZW5kCiMjIyMgU3RvbHR6Y29pbiBmbHVjdHVhdGVzIHRocm91Z2hvdXQgdGhlIHllYXIKCmBgYHtyfQpnb2luZ191cF9ieSA9IDIKdHJlbmQgPSBzZXEoZnJvbSA9IDAsIAogICAgICAgICAgICB0byA9IChuLTEpICogZ29pbmdfdXBfYnksCiAgICAgICAgICAgIGJ5ID0gZ29pbmdfdXBfYnkpCnNlYXNvbmFsID0gc2luKDE6bikgKiAxMAoKdHNfbGluX3RyZW5kX3NlYXNvbmFsID0gdGliYmxlKGRhdGUgPSBkYXRlcywgc3RvbHR6Y29pbiA9IHRpbWVfc2VyaWVzJHN0b2x0emNvaW4gKyB0cmVuZCArIHNlYXNvbmFsLCAKICAgICAgICAgICAgICAgICAgICAgdHlwZSA9ICdub24tc3RhdGlvbmFyeScsIHRyZW5kID0gJ2xpbmVhcicpCgpnZ3Bsb3QodHNfbGluX3RyZW5kX3NlYXNvbmFsLCBhZXMoeCA9IGRhdGUsIHkgPSBzdG9sdHpjb2luKSkgKyAKICBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgc2UgPSBGQUxTRSkgKyAKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gbWVhbihzdG9sdHpjb2luKSksIGNvbG9yPSJyZWQiKQpgYGAKCgoKCiMjIyMgUGVyZm9ybSByZWdyZXNzaW9uIGFuZCBsb29rIGF0IHRoZSByZXNpZHVhbHMKCmBgYHtyfQp0c19maXQgPSBsbShzdG9sdHpjb2luIH4gZGF0ZSwgZGF0YSA9IHRzX2xpbl90cmVuZF9zZWFzb25hbCkKdHNfZml0CmBgYAoKCmBgYHtyfQpyZXMgPSB0aWJibGUocmVzaWR1YWxzID0gcmVzaWQodHNfZml0KSkKZ2dwbG90KHJlcywgYWVzKHggPSAxOm5yb3cocmVzKSwgeSA9IHJlc2lkdWFscykpICsgZ2VvbV9wb2ludCgpICsgZ2VvbV9ydWcoc2lkZXMgPSAnbCcpCmBgYAoKYGBge3J9CmdncGxvdChyZXMsIGFlcyh4ID0gcmVzaWR1YWxzKSkgKyBnZW9tX2hpc3RvZ3JhbShhZXMoeSA9IC4uZGVuc2l0eS4uKSwgYmlucyA9IDQwKSArIGdlb21fZGVuc2l0eShjb2wgPSAncmVkJywgc2l6ZSA9IDEuNSkKYGBgCgoKYGBge3J9CmFkZi50ZXN0KHpvbzo6em9vKHRzX2xpbl90cmVuZF9zZWFzb25hbCRzdG9sdHpjb2luLCB0c19saW5fdHJlbmRfc2Vhc29uYWwkZGF0ZSkpCmBgYAoKCgoKIyMjIEFkZCBhbiBleHBvbmVudGlhbCB0cmVuZCB0byB0aGUgc2Vhc29uYWwgY29tcG9uZW50CiMjIyMgU3RvbHR6Y29pbiBoYXMgaW5jcmVhc2VkIGluIHZvbGF0aWxpdHkgb3ZlciB0aGUgeWVhcnMKICAqIFN0b2x0emNvaW4gbm93IGhhczoKICAgICogQSBwb3NpdGl2ZSBhbmQgZ3Jvd2luZyBsaW5lYXIgdHJlbmQKICAgICogU2Vhc29uYWxpdHkKCmBgYHtyfQpnb2luZ191cF9ieSA9IDIKdHJlbmQgPSBzZXEoZnJvbSA9IDAsIAogICAgICAgICAgICB0byA9IChuLTEpICogZ29pbmdfdXBfYnksCiAgICAgICAgICAgIGJ5ID0gZ29pbmdfdXBfYnkpCnNlYXNvbmFsID0gc2luKDE6bikgKiAxOm4KCnRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAgPSB0aWJibGUoZGF0ZSA9IGRhdGVzLCBzdG9sdHpjb2luID0gdGltZV9zZXJpZXMkc3RvbHR6Y29pbiArIHRyZW5kICsgc2Vhc29uYWwsIAogICAgICAgICAgICAgICAgICAgICB0eXBlID0gJ25vbi1zdGF0aW9uYXJ5JywgdHJlbmQgPSAnbGluZWFyJykKCmdncGxvdCh0c19saW5fdHJlbmRfc2Vhc29uYWxfZXhwLCBhZXMoeCA9IGRhdGUsIHkgPSBzdG9sdHpjb2luKSkgKyAKICBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCAKICAgICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgc2UgPSBGQUxTRSkgKyAKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gbWVhbihzdG9sdHpjb2luKSksIGNvbG9yPSJyZWQiKQpgYGAKCgoKCiMjIyMgUGVyZm9ybSByZWdyZXNzaW9uIGFuZCBsb29rIGF0IHRoZSByZXNpZHVhbHMKCmBgYHtyfQp0c19maXQgPSBsbShzdG9sdHpjb2luIH4gZGF0ZSwgZGF0YSA9IHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHApCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCmBgYHtyfQphZGYudGVzdCh6b286Onpvbyh0c19saW5fdHJlbmRfc2Vhc29uYWxfZXhwJHN0b2x0emNvaW4sIHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpCmBgYAoKCiMjIFIncyBkZWNvbXBvc2UgZnVuY3Rpb24gZG9lcyBhIGxvdCBmb3IgeW91CmBgYHtyfQpkMXkgPSAgYXMubnVtZXJpYyhsdWJyaWRhdGU6OnllYXIobWluKHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQpkMW0gPSBhcy5udW1lcmljKGx1YnJpZGF0ZTo6bW9udGgobWluKHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQpkMnkgPSAgYXMubnVtZXJpYyhsdWJyaWRhdGU6OnllYXIobWF4KHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQpkMm0gPSBhcy5udW1lcmljKGx1YnJpZGF0ZTo6bW9udGgobWF4KHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkZGF0ZSkpKQoKeCA9IHRzKHRzX2xpbl90cmVuZF9zZWFzb25hbF9leHAkc3RvbHR6Y29pbiwgc3RhcnQ9YyhkMXksIGQxbSksIGVuZD1jKGQyeSwgZDJtKSwgZnJlcXVlbmN5PTEyKQoKcGxvdChkZWNvbXBvc2UoeCwgdHlwZSA9ICdhZGRpdGl2ZScpKQoKYGBgCgoKYGBge3J9CnBsb3QoZGVjb21wb3NlKHgsIHR5cGUgPSAnbXVsdGlwbGljYXRpdmUnKSkKYGBgCgoKIyMjIEV4YW1wbGUgZnJvbSBkZWNvbXBvc2UgZnVuY3Rpb24KIyMjIyByYW5kb20gcG9ydGlvbiBzaG91bGQganVzdCBsb29rIGxpa2Ugd2hpdGUgbm9pc2UgKG90aGVyd2lzZSB5b3UndmUgbGlrZWx5IG1pc3NlZCBzb21ldGhpbmcpCmBgYHtyfQpyZXF1aXJlKGdyYXBoaWNzKQoKbSA8LSBkZWNvbXBvc2UoY28yKQpwbG90KG0pCmBgYAoKCgpgYGB7ciwgd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRX0KZG93ID0gdHFfaW5kZXgoIkRPVyIpICU+JQogIHRxX2dldChnZXQgPSAic3RvY2sucHJpY2VzIikKaGVhZChkb3csIDEwKQpgYGAKCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFfQpwID0gZ2dwbG90KGRvdywgYWVzKHggPSBkYXRlLCB5ID0gY2xvc2UsIGNvbCA9IHN5bWJvbCkpICsgCiAgZ2VvbV9saW5lKCkgKyAKICBnZ2hpZ2hsaWdodChtYXgoY2xvc2UpID4gMjUwKSArIAogIHRoZW1lX21pbmltYWwoKSArIAogIGZhY2V0X3dyYXAofiBzeW1ib2wpCnByaW50KHApCmBgYAoKCgojIyBJcyBCQSBhIHN0YXRpb25hcnkgdGltZSBzZXJpZXM/CgojIyMjIENhbiB5b3Ugc2ltcGx5IHpvb20gaW4gYW5kIGZpbmQgYSBwb3J0aW9uIHRoYXQgaXMgc3RhdGlvbmFyeT8KYGBge3J9CmJhID0gZG93ICU+JSAKICBmaWx0ZXIoc3ltYm9sID09ICdCQScpICU+JSAKICBzZWxlY3QoZGF0ZSwgY2xvc2UpCnBsb3RseTo6Z2dwbG90bHkoZ2dwbG90KGJhLCBhZXMoeCA9IGRhdGUsIHkgPSBjbG9zZSkpICsgZ2VvbV9saW5lKCkgKyBnZW9tX3Ntb290aChtZXRob2QgPSBsbSwgc2UgPSBGQUxTRSkpCmBgYAoKCmBgYHtyfQphZGYudGVzdCh6b286OnpvbyhiYSRjbG9zZSwgYmEkZGF0ZSkpCmBgYAoKCgpgYGB7cn0KdHNfZml0ID0gbG0oY2xvc2UgfiBkYXRlLCBkYXRhID0gYmEpCnRzX2ZpdApgYGAKCgpgYGB7cn0KcmVzID0gdGliYmxlKHJlc2lkdWFscyA9IHJlc2lkKHRzX2ZpdCkpCmdncGxvdChyZXMsIGFlcyh4ID0gMTpucm93KHJlcyksIHkgPSByZXNpZHVhbHMpKSArIGdlb21fcG9pbnQoKSArIGdlb21fcnVnKHNpZGVzID0gJ2wnKQpgYGAKCmBgYHtyfQpnZ3Bsb3QocmVzLCBhZXMoeCA9IHJlc2lkdWFscykpICsgZ2VvbV9oaXN0b2dyYW0oYWVzKHkgPSAuLmRlbnNpdHkuLiksIGJpbnMgPSA0MCkgKyBnZW9tX2RlbnNpdHkoY29sID0gJ3JlZCcsIHNpemUgPSAxLjUpCmBgYAoKCgoKCgojIyMjIEFyZSB0aGUgZGFpbHkgcmV0dXJucyBzdGF0aW9uYXJ5PwoKYGBge3J9CnJldCA9IGJhICU+JQogIHRxX3RyYW5zbXV0ZShjbG9zZSwgbXV0YXRlX2Z1biA9IGRhaWx5UmV0dXJuKQpnZ3Bsb3QocmV0LCBhZXMoeCA9IGRhdGUsIHkgPSBkYWlseS5yZXR1cm5zKSkgKyBnZW9tX2xpbmUoKSArIGdlb21fc21vb3RoKG1ldGhvZCA9IGxtLCBzZSA9IEZBTFNFKQpgYGAKCgpgYGB7cn0KYWRmLnRlc3Qoem9vOjp6b28ocmV0JGRhaWx5LnJldHVybnMsIHJldCRkYXRlKSkKYGBgCgoKYGBge3J9CnRzX2ZpdCA9IGxtKGRhaWx5LnJldHVybnMgfiBkYXRlLCBkYXRhID0gcmV0KQp0c19maXQKYGBgCgoKYGBge3J9CnJlcyA9IHRpYmJsZShyZXNpZHVhbHMgPSByZXNpZCh0c19maXQpKQpnZ3Bsb3QocmVzLCBhZXMoeCA9IDE6bnJvdyhyZXMpLCB5ID0gcmVzaWR1YWxzKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX3J1ZyhzaWRlcyA9ICdsJykKYGBgCgpgYGB7cn0KZ2dwbG90KHJlcywgYWVzKHggPSByZXNpZHVhbHMpKSArIGdlb21faGlzdG9ncmFtKGFlcyh5ID0gLi5kZW5zaXR5Li4pLCBiaW5zID0gNDApICsgZ2VvbV9kZW5zaXR5KGNvbCA9ICdyZWQnLCBzaXplID0gMS41KQpgYGAKCgpgYGB7cn0KCmBgYAoK